home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / four.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  7.5 KB  |  301 lines
  1. /* standard begin
  2. */
  3.  
  4. #include <stdio.h>
  5. #include <spec.h>
  6.  
  7. #define  MAXPEAKS 10
  8.  
  9. float x,y,*spc,*err,*tim,*tr,*ti;
  10.  
  11. help()
  12. {
  13.   printf("Fast fourier transformation from input spectrum\n");
  14.   printf("The resulting spectrum is allready a power spectrum\n");
  15.   printf("option     -iv        inverse spectrum\n");
  16.   printf("           -rf name   save real part to file\n");
  17.   printf("           -if name   save imaginary part to file\n");
  18.   printf("           -pc n      append n peaks to comment\n");
  19.   printf("           -pp n      print n peaks to stdout\n");
  20.   printf("           -s         suppress printing of powerspectrum to stdout\n");
  21. }
  22.  
  23. main(argc,argv)
  24. int argc;
  25. char *argv[];
  26. {
  27. int kmax2,n,max,inv,peaks[MAXPEAKS];
  28. char z[1024],comment[1024];
  29.  
  30.    inv=0;
  31.    if(checkopt(argc,argv,"-iv",z)) inv=1;
  32.  
  33.    tr = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  34.    ti = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  35.    if(ti==NULL) {
  36.       printf("sorry, not enough memory\n");
  37.       exit(-1);
  38.    }
  39.    for(n=0;n<=_MAXSPCLEN;n++) {  /* should not be neccessary, but good style */
  40.       tr[n]=0; ti[n]=0;
  41.    }
  42.    spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  43.    err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  44.    tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
  45.    if(tim==NULL) {
  46.       printf("sorry, not enough memory\n");
  47.       exit(-1);
  48.    }
  49.  
  50.    max=readspec(argv[1],spc,err,tim,comment);
  51.    for(n=0;n<=max;n++) tr[n]=spc[n];
  52.    kmax2=four(inv,max);
  53.    for(n=0;n<=max;n++) {         /* now recalculate time spectrum */
  54.       x = (float) n;
  55.       tim[n] = (6.2831853 * ((float) n)) / (((float) kmax2) * _tica);
  56.    }
  57.    for(n=0;n<=max;n++) spc[n]=sqrt((tr[n]*tr[n])+(ti[n]*ti[n]));
  58.    fpeaks(peaks,MAXPEAKS,spc,kmax2/2);
  59.    if(checkopt(argc,argv,"-pc",z)) {
  60.       n=atoi(z);
  61.       peaktoa(z,n,peaks,kmax2,spc);
  62.       strcat(comment,": ");
  63.       strcat(comment,z);
  64.    }
  65.    if(!checkopt(argc,argv,"-s",z)) writespec("",spc,err,kmax2/2,2,comment);
  66.    if(checkopt(argc,argv,"-rf",z)) writespec(z,tr,err,kmax2/2,2,comment);
  67.    if(checkopt(argc,argv,"-if",z)) writespec(z,ti,err,kmax2/2,2,comment);
  68.    if(checkopt(argc,argv,"-pp",z)) {
  69.       n=atoi(z);
  70.       peaktoa(z,n,peaks,kmax2,spc);
  71.       printf("%s\n",z);
  72.    }
  73.    if(checkopt(argc,argv,"-o",z)) {            /* write modified time array */
  74.       strcat(z,".tim");
  75.       writespec(z,tim,err,kmax2/2,2,comment);
  76.    }
  77.    free(tr); free(ti); free(spc); free(err); free(tim);
  78.    exit(0);
  79. }
  80.  
  81. /* ------------------------------------------------------------------------
  82.    Fast Fourier transform
  83.    ------------------------------------------------------------------------ */
  84. log2(n)
  85. int n;
  86. {
  87. int erg;
  88.  
  89.    erg=0;
  90.    while((1<<erg)<n) erg++;
  91.    return(erg);
  92. }
  93.  
  94. four(inv,max)
  95. int inv;
  96. {
  97. int   k,ii,i,j,it,i0,i1,i2,i3,
  98.       kk,j0,kk1,l,l1,jp1,kp1,j1,j2,
  99.       kmax2;
  100. float z,zr,zi,wr,wi,ws,ur[15],ui[15],um,th;
  101.  
  102.    it=max;
  103.    um=0.5;
  104.    for(i=1;i<=15;i++) {
  105.       um = 0.5 * um;
  106.       th = 6.283185 * um;
  107.       ur[i] = cos(th); ui[i]=sin(th);
  108.    }
  109.    um = 1 - (2*inv);
  110.    i0 = log2(it); kmax2 = 1 << i0;
  111.    it = kmax2;
  112.    ii=i0 ; i1=it/2; i3=1;
  113.    do {
  114.       k=0; i2=i1+i1;
  115.       do {
  116.          wr=1; wi=0; kk=k; j0=i0;
  117.          while(kk!=0) {
  118.             do {
  119.                j0=j0-1; kk1=kk; kk=kk/2;
  120.             } while((kk1-(2*kk)) == 0);
  121.             ws=(wr*ur[j0]) - (wi*ui[j0]);
  122.             wi=(wr*ui[j0]) + (wi*ur[j0]);
  123.             wr = ws;
  124.          }
  125.          wi=wi*um; j=0;
  126.          do {
  127.             l=j*i2+k+1; l1=l+i1;
  128.             zr=tr[l]+tr[l1];
  129.             zi=ti[l]+ti[l1];
  130.             z=wr*(tr[l]-tr[l1]) - (wi*(ti[l]-ti[l1]));
  131.             ti[l1]=wr*(ti[l]-ti[l1]) + (wi*(tr[l]-tr[l1]));
  132.             tr[l]=zr; tr[l1]=z; ti[l]=zi; j=j+1;
  133.          } while((j-i3)<0);
  134.          k=k+1;
  135.       } while((k-i1)<0);
  136.       i3 = i3 + i3; i0 = i0 - 1; i1 = i1 / 2;
  137.    } while(i1 > 0);
  138.    j=1; um=1;
  139.    if(inv==1) um=1/it;
  140.    do {
  141.       k=0; j1=j;
  142.       for(i=1;i<=ii;i++) {
  143.          j2 = j1 / 2; k=(2*(k-j2))+j1; j1=j2;
  144.       }
  145.       jp1=j+1;
  146.       if((k-j) >= 0) {
  147.          if((k-j) == 0) {
  148.             tr[jp1] = um * tr[jp1];
  149.             ti[jp1] = um * ti[jp1];
  150.          } else {
  151.             kp1 = k + 1;
  152.             zr = tr[jp1]; zi=ti[jp1];
  153.             tr[jp1] = um * tr[kp1];
  154.             ti[jp1] = um * ti[kp1];
  155.             tr[kp1] = zr * um;
  156.             ti[kp1] = zi * um;
  157.          }
  158.       }
  159.       j = jp1;
  160.    } while((j-it+1)<0);
  161.    tr[1] = um * tr[1];
  162.    ti[1] = um * ti[1];
  163.    return(kmax2);
  164. }
  165.  
  166. /* -----------------------------------------------------------------------
  167.    find peaks in powerspectrum
  168.    ----------------------------------------------------------------------- */
  169.  
  170. fpeaks(peaks,pmax,spc,smax)
  171. float spc[];
  172. int peaks[],pmax,smax;
  173. {
  174. int n,sptr;
  175.  
  176.    sptr=0;
  177.    for(n=0;n<=pmax;n++) peaks[n]= -1 ;             /* mark peaks as unset */
  178.    for(;;) {
  179.       while(spc[sptr] <= spc[sptr+1]) {             /* find next maximum */
  180.          sptr=sptr+1;
  181.          if(sptr>smax) return(0);
  182.       }
  183.       peaks[pmax]=sptr;                            /* smallest peak redefined */
  184.       peaksort(peaks,pmax,spc,smax);               /* sort peaks for height  */
  185.       while(spc[sptr] > spc[sptr+1]) {             /* find next minimum */
  186.          sptr=sptr+1;
  187.          if(sptr>smax) return(0);
  188.       }
  189.    }
  190. }
  191.  
  192. peaksort(peaks,pmax,spc,smax)
  193. float spc[];
  194. int peaks[],pmax,smax;
  195. {
  196. int m,p1,p2,flag;
  197.  
  198.    flag=0;
  199.    if(peaks[pmax]==-1) return(0);
  200.    while(flag == 0) {
  201.       flag=1;
  202.       for(m=0;m<pmax;m++) {
  203.          p1=peaks[m];
  204.          p2=peaks[m+1];
  205.          if(p1 == -1) {
  206.             peaks[m]=pmax;
  207.         flag=0;
  208.             continue;
  209.          } 
  210.          if(p2 == -1) {
  211.         peaks[m+1]=pmax;
  212.         flag=0;
  213.         continue;
  214.          }
  215.          if(spc[p1]<spc[p2]) {
  216.             peaks[m]=p2;
  217.             peaks[m+1]=p1;
  218.         flag=0;
  219.          }
  220.       }
  221.    }
  222. }
  223.  
  224. /* -------------------------------------------------------------- */
  225.  
  226. float peakmean(peak,spc)
  227. int peak;
  228. float spc[];
  229. {
  230. int n,left,right;
  231. float x,y,y1,y2,sum,lsum,mean;
  232.  
  233.    left=peak;                              /* finding left side of peak */
  234.    for(n=0;n<=20;n++) {
  235.       if(left==0) break;
  236.       y1=spc[left];
  237.       y2=spc[left-1];
  238.       if(y2>y1) break;
  239.       y = y1 / 10.0;
  240.       x = fabs(y1-y2);
  241.       if(x<y) break;
  242.       left = left - 1;
  243.    }
  244.    right=peak;                             /* finding right side of peak */
  245.    for(n=0;n<=20;n++) {
  246.       y1=spc[right];
  247.       y2=spc[right+1];
  248.       if(y2>y1) break;
  249.       y = y1 / 10.0;
  250.       x = fabs(y1-y2);
  251.       if(x<y) break;
  252.       right = right + 1;
  253.    }
  254.    sum = 0.0;                                /* peak integral */
  255.    for(n=left;n<=right;n++) sum = sum + spc[n];
  256.    sum = sum / 2;
  257.    lsum = 0.0;
  258.    n=left;
  259.    while(lsum < sum) {                       /* finding half integral */
  260.       y1 = lsum;
  261.       lsum = lsum + spc[n++];
  262.    }
  263.    n = n - 1;
  264.    y1 = y1 - sum;
  265.    y2 = lsum - sum;
  266.    if(y1==y2) {
  267.       x = peak;
  268.       return(x);
  269.    }
  270.    x = y1 / (y1 - y2);
  271.    mean = (n-1) + x;
  272.    return(mean);
  273. }
  274.  
  275. peaktoa(s,max,peaks,kmax2,spc)
  276. char s[];
  277. int max,peaks[],kmax2;
  278. float spc[];
  279. {
  280. char z[80];
  281. int i,n,m,ix;
  282. float x,x1,x2;
  283.  
  284.    strcpy(s,"");
  285.    for(n=0;n<max;n++) {
  286.       i=peaks[n];                                 /* get maximum */
  287.       x=peakmean(i,spc);                          /* mean instead of maximum */
  288.       ix=x;
  289.       x1=tim[ix]; x2=tim[ix+1];                   /* interpolation between    */
  290.       x = x1 + ((x-ix) * (x2-x1));                /* tim[i] and tim[i+1]      */
  291.       sprintf(z," %E",x);
  292.       if(fabs(x)>0.001) sprintf(z," %-1.5f",x);
  293.       if(fabs(x)>1.0) sprintf(z," %-3.2f",x);
  294.       for(m=0;m<n;m++) {
  295.      if(peaks[m] == i) return(0);
  296.       }
  297.       strcat(s,z);
  298.    }
  299. }
  300.  
  301.